Liquid crystal director fluctuations and surface anchoring by molecular simulation 
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We propose a simple and reliable method to measure the liquid crystal surface anchoring strength 
by molecular simulation. The method is based on the measurement of the long-range fluctuation 
modes of the director in confined geometry. As an example, molecular simulations of a liquid crystal 
in slab geometry between parallel walls with homeotropic anchoring have been carried out using the 
Monte Carlo technique. By studying different slab thicknesses, we are able to calculate separately 
the position of the elastic boundary condition, and the extrapolation length. 
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I. INTRODUCTION 

Liquid crystal anchoring effects have been intensively 
studied both experimentally and theoretically during the 
past decades jl], |). Such an interest is well understood 
since most liquid crystal devices are cells comprising ori- 
enting surfaces with a liquid crystal in between. Typi- 
cally, aligning surfaces provide a uniform orientation of 
the liquid crystal director in the cell interior. 

On the phenomenological level, liquid crystal anchor- 
ing can be described by two basic parameters: the easy 
axis direction, e, and the anchoring coefficient W . These 
two parameters are critical design parameters for every 
liquid crystal device || . A variety of experimental meth- 
ods have been proposed to measure anchoring param- 
eters, in particular the anchoring coefficient W [|[ |). 
Most of them measure surface director deviations in an 
external field and involve rather complicated optical se- 
tups. 

In spite of the practical importance, the mechanism 
of the director alignment is still not well understood. 
Experimental techniques always involve optical measure- 
ments. They test the entire liquid crystal cell and there- 
fore cannot provide a satisfactory description of the thin 
interface region. Theoretical investigations of anchoring 
are also quite controversial. For example, the usual con- 
tinuous phenomenological theory has divergent surface 
terms in the elastic free energy expansion ||, |7j. 

One of the approaches for the systematic investigation 
of anchoring phenomena is computer simulation of liquid 
crystals in confined geometries. Indeed, computer simu- 
lation is a well established method to treat bulk elastic 
coefficients |j, [| , the surface anchoring strength , and 
structures of disclination cores |12|| . This means that 
computer simulation allows investigation of details of the 
liquid crystalline structure which cannot be resolved ex- 
perimentally 

Several papers have already been published trying to 
search for reasonable surface potentials to use in simula- 
tions (H, H, H |||. However, most of them do 
not characterize aligning surfaces using well established 
parameters. Questions about the formation of a solid in- 



terface layer, values of the easy axis angle and anchoring 
coefficient are still open. The reason for this is probably 
a lack of reliable methods to measure these parameters 
by computer simulation. 

In this paper, we propose a technique to measure the 
surface anchoring strength W by computer simulation. 
The method itself is based on the study of the director 
fluctuations in the liquid crystal cell. A similar approach 
has already been used for the experimental characteri- 
zation of the interface |]l9|| and is based on measuring 
the light scattering caused by the director fluctuations. 
In computer simulation, the director fluctuations can be 
studied directly, using ensemble averages of functions of 
the second-rank order tensor components. A fit to the 
fluctuation amplitudes with equations predicted by elas- 
tic theory then allows determination of the surface an- 
choring strength. 



II. THEORY 

Large length- and time-scale fluctuations of the direc- 
tor n(r) can be described in the continuum model of liq- 
uid crystals, based on the phenomenological elastic con- 
stants. In this approach, the hydrodynamic equations for 
the director and the boundary conditions can be obtained 
by minimization of the cell free energy EfJ : 



F = F h + F h + F s 



(1) 



where Fy, — ^j(dn/dt) 2 is a hydrodynamic term with 
an effective viscosity coefficient 7, Ft, is the Frank elastic 
free energy, and F s is the surface anchoring energy. 

In what follows we use the one-elastic-constant approx- 
imation, i.e. K\\ = K22 = K33 = K. Then the Frank 
free energy can be brought into the form: 



F h = Ik 



(V • n) 2 + (V x n) 



dr . 



The liquid crystal cell is bounded by surfaces z = 0, L 
which provide some kind of anchoring condition jl| . Be- 
low we consider homeotropic anchoring, that is, normal 
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to the surface. Planar anchoring can be treated in the 
same way. We assume that the interaction of the director 
with the cell surfaces has the Rapini-Papoular form JnJ : 



F s = --W 



SQ;Sl 



(n-e ,i) dr ± , 



where unit vectors e £ define directions of the easy axes 
at z = 0,L: e = e 2 , = — e 2 ; W is the anchoring 
energy, and rj_ = (x,y). 

Equations for the director and boundary conditions 
can be obtained by minimization of the cell free energy 
eqn (|l|). In the case of homeotropic anchoring on both 
surfaces the stationary director distribution in the cell 
is homeotropic, i.e. no = e 2 . Therefore, we have to in- 
vestigate small perturbations of the director around the 
distribution: 



n(r) = no + <5n(r) 



(2) 



Minimizing the total free energy (Q) and linearizing the 
equations for the director and boundary conditions with 
respect to <5n, we obtain equations 



and boundary conditions 



K V 2 8n . 



d 

WSn±K—5n 

oz 



(3) 



(4) 



-L.Q 



for the fluctuations. Taking into account eqn. (g), the ex- 
pression for the free energy fluctuations can be rewritten 
in the form of the average of the self-conjugate (Hermi- 
tian) operator — ^KV 2 : 



SF h = ~K [ 5n(r)V 2 5n(r) dr . 
2 Jv 



(5) 



The eigenfunctions of the operator — iifV 2 , which sat- 
isfy boundary conditions (jij), form a complete set of 
orthogonal functions characterized by wave vector q. 
Therefore, <5n(r) can be expanded in a series of these 
orthogonal functions: 

Sn(r) = i ^ e icur± x 



where = (q x ,q y ), and 



(T) 



Here we have introduced the dimensionless wave vector 
X = q z L and anchoring parameter 



WL 



L 

A ' 



(8) 



where A is the extrapolation length (g0|. The wave vec- 
tors q z form a discrete spectrum because of confinement 
in the z direction which depends on the anchoring energy 
W (see Appendix for details). The explicit form of the 
q z spectrum is given by the secular equation: 



(£ 2 -X 2 ) sinx + 2£xcosx = 0. 



(9) 



Each individual mode can now be seen from eqn Q3J) to 
relax exponentially with a relaxation time given by 



j/K (q 2 ± + q 2 z ) 



(10) 



Substituting expansion into the free energy (|J) , and 
performing the integration over the cell volume we ob- 
tain: 



SF b 



K y, q 2 (2t + X 2 +e) x 



V 



(ix + C) 2 



W+)(qx,<Z z )W + )(-qj 



(11) 



Integrating, we took into account the orthogonality of the 
eigenfunctions in the expansion (^|) with different eigen- 
vectors q, which allowed us to reduce the summations 
over q, q' to a single sum over q. 

Application of the equipartition theorem of classical 
statistical mechanics, just as for elastic fluctuations in 
bulk liquid crystals p2], gives the fluctuation amplitudes: 



W+>(q ± ,g 2 )<5nM(-q ± ,^)) = 

kBTV(i X + Q 2 
2Kq 2 (2£ + x 2 + e) ' 



(12) 



where (...) denotes an ensemble average. 

In molecular simulations, rather than measuring direc- 
tor fluctuations, it is more convenient to measure fluctua- 
tions of the second-rank order tensor components (follow- 
ing Forster [p2[ ). We define the real-space order tensor 
density 

V \ - 

Q a p(r) = — 2^*(r-ri)Q^, 



3 ( 1 x 

la/3 = 2 I UiaUi P ~ Y a P 



where a,/3 = x, y, z, in terms of the orientation vectors 
of each molecule i (we consider only uniaxial molecules). 
If we assume that there is no variation in the degree of 
ordering, we may write 

3 1 

Qap(r) = ^Qn a (r)n (r) - -Q5 af3 

where Q is the order parameter, i.e. the largest eigenvalue 
of Qa/3 (r) ■ If the director n is taken to lie along the z 
axis throughout the sample, the off-diagonal components 
*9o!z(r), a = x,y, are proportional to the fluctuations of 
the corresponding director components 
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This is the situation for homeotropic anchoring at both x and the surface normal along z, the components Q xy , 
surfaces. For planar anchoring, with the director along Q xz are important (and non-equivalent). 

I 



Measurements are performed directly in reciprocal space. The Fourier transform of the real-space order tensor is 

Q Q/J (k) = jf Q a/3 (r)e ik ' r dr = ^ ]T Q^e ikr * . 



Then the fluctuations ( |Qd/3 (k) | 2 \ can be easily measured from simulations: 



q\p cos(k • r< ) ) qU sin ( k ■ r *) 



YL 

N 2 



(13) 



We explicitly relate simulation-measured fluctuation modes with theoretically predicted amplitudes of the director 
fluctuations for = 

Q az (k z ) = Sn^(0,q z ) 

where k = k z L. Note that the q z take discrete (but not equally spaced) values as discussed earlier, while the k z values 
are unrestricted. 

Since i5n(r) is real, using eqn ([|) we have 







e i(«"X) _ 1" 


K + X 


W+u 


K-X 



e+x 2 



■6n^(-q ± ,q z ) 



Taking this equation into account, and the fact that fluctuations with different wavevectors are independent, i.e. 
(Jn( + ) (q_L, q z )Sn^ (— q^, q' z )) = if q z =/= q' z , the corresponding ensemble average of the squared order parameter 
can be rewritten as 



Q az (k z )\ 2 ) = %T^-Yl 



x 2 +e 



k ^ q 2 z m+x 2 +e) 



K + X 



£\ e K«-x) _ i 



(14) 



We measure Q and \ \Q a z(k z )\ ) from simulations, eqn (|13|), and then compare with the theoretical prediction, 
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eqn (14), which is parametrized by L, A and K. Both the permitted q z spectrum, and the variation of ^|Q QZ (fc 

with k z are sensitive to the anchoring strength parameter £ = L/ A . 

Fluctuation amplitudes given by eq. ( |l4] ) have features that simplify the fitting procedure. First, terms with small 
q z values dominate because of the q z in the denominator of eq (jl4|). Therefore, it is always possible to truncate this 
sum and use only the first values of the q z spectrum. Then, for large k z or for k >> x, eq. ( |l4| ) can be simplified so 
that the dependence on k z is explicit: 



\QM\ 2 ) = \^T^- 



sin (k/2) 
k/2 



E 



X 



q 2 W + X 2 + e)' 



(15) 



Therefore, for large k z , the fluctuation amplitude {^Q a z{k z )\ ) has a characteristic oscillation with the period 
given by k = k z L = 2ir. This means that we can adjust the cell thickness L independently of parameters A and K by 
examining characteristic wavelength of the fluctuation amplitude \\Q a z(k z )\ 2 

I 



III. MOLECULAR MODEL AND SIMULATION 
METHODS 

To test the technique proposed, we simulated a liquid 
crystal confined between parallel walls (slab geometry), 



with finite homeotropic anchoring at the walls. The di- 
rector fluctuations occur around the preferred, uniform, 
alignment perpendicular to the walls. 
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We performed Monte Carlo (MC) simulation of the liq- 
uid crystal system. We used a molecular model which has 
been studied earlier in this geometry jl0| . The molecules 
in this study were modelled as hard ellipsoids of revolu- 
tion of elongation e = a/6 = 15, where a is the length 
of the semi-major axis and 6 the length of the two equal 
semi-minor axes. The phase diagram and properties of 
this family of models are well studied @ |], H EH 13. 
It is useful to express the density as a fraction of the 
close-packed density p cp of perfectly aligned hard ellip- 
soids, assuming an affinely-transformed face-centred cu- 
bic or hexagonal close packed lattice. In this case, the 
isotropic-nematic phase transition occurs at quite a low 
density, p/p cp ~ 0.2, and the simulations are performed 
at a state point corresponding to pj p cp ~ 0.28, for which 
the nematic order parameter is Q ~ 0.85. For this model, 
temperature is not a significant thermodynamic quantity, 
so it is possible to choose k^T = 1 throughout. 

The slab geometry is defined by two hard parallel con- 
fining walls, which cannot be penetrated by the centres 
of the ellipsoidal molecules. Packing considerations gen- 
erate homeotropic ordering at the surface. Surface an- 
choring in this system has been studied recently |l(| for 
a system with wall separation L z — 1256 = 8.33a, by 
applying an orienting perturbation at one of the walls 
and observing the response at the other. This yielded an 
estimate of the extrapolation length A ~ 356 w 2.33a. 
In the current work, simulations were carried out for 
systems of N = 2000 particles with wall separations 
L z = 6.58a, 8.22a, 9.86a, which (from the above estimate 
of A) would correspond to surface anchoring parameters 
in the range 2.8 < £ < 4.2. 

To be sure that we have the same state point for each 
simulation, we adjusted the density to have the same P zz 
component of the pressure tensor for all systems. Then 
a sequence of runs was carried out using the constant- 
ly! 1 ensemble, allowing typically 10 5 MC sweeps for 
equilibration and 10 7 sweeps for accumulation of averages 
(one sweep is one attempted move per particle). 



IV. SIMULATION RESULTS AND DISCUSSION 

The simulation results were analysed to give a density 
profile, and an order tensor profile, which are shown in 
Figs |l|, |2| From these profiles we can see that the walls 
are sufficiently well separated, and the variation of the 
order parameter across the slab is small, even in spite of 
the large change in local density near the walls. 

The order tensor fluctuations in reciprocal space were 
calculated using expression (|l3|). To fit the simulation 
results with the elastic theory we have to remember that 
the size of the simulation box L z is not necessarily equal 
to the liquid crystal cell thickness L appearing in the 
elastic theory. The former is a physical quantity, which, 
in statistical mechanics, is determined by the positions 
at which the liquid number density becomes identically 
zero. The latter appears in the elastic theory: it is deter- 
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FIG. 1: Density profiles as functions of z-coordinate. Dis- 
tances are normalized by the molecular semi-axis length a and 
the density is expressed relative to the closed-packed density 
p cp for perfectly aligned ellipsoids. The profiles are symmetri- 
cal, only one side of the slab is shown. The results for different 
wall separations are almost indistinguishable on this scale. 

0.9 1 1 1 1 1 1 1 1 1 1 1 
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FIG. 2: Profiles of nematic order parameter Q for different 
wall separations. Distances are normalized by the molecular 
semi-axis length a. The profiles are symmetrical, only one 
side of the slab is shown. Note the highly expanded vertical 
scale. 



mined by the positions at which the orientational elastic 
boundary conditions are applied. Physically, the differ- 
ence may be ascribed to partial penetration of the walls 
by the liquid crystal molecules, formation of an ordered 
(or solid) layer near the walls, or other molecular-scale 
features. We assume that we may write L = L z + 2L W , 
where the value L w (which may be positive or negative) 
is characteristic of the wall, independent of L z , and may 
be determined in our fitting process. 

The best estimate of the wall-induced separa- 
tion distance L w was obtained examining the ra- 
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FIG. 3: Ratio of the director fluctuations as function of 
wavevector (normalized by the molecular semi-axis length a) 
for different wall separations. Symbols: Monte Carlo results; 
solid lines: elastic theory; dashed lines: elastic theory, without 
correction for the difference between the elastic-theory cell 
thickness L and the simulation box size L z . 



tios ^|Q QZ (fc z ,ii)| 2 ^> / (\Q az (k Zl L 2 )\ 2 ^ for different L lj2 

since they reveal more structure for large k z . Plotting the 
results in this way removes the overall scaling of the am- 
plitudes of fluctuations, which are sensitive to changes in 
A, while the shapes of the curves, and the characteristic 
oscillation wavelengths, are sensitive to the choice of i w . 
These ratios with L\ = 6.58a, 8.22a, L 2 = 9.86a and cor- 
responding fitting curves are plotted in Fig. |. The best 
fits were obtained with L w /a = 0.59. In this figure we 
also plot theoretical curves with L w = 0: the clear dis- 
crepancy with the simulation results indicates that the 
simulation box size L z is indeed significantly different 
from the actual cell thickness. 

Following the determination of L w , we adjusted the 
extrapolation length A to obtain the best fit to the fluc- 
tuation data: the fluctuation amplitudes with small k z 
are most sensitive to this quantity. Together with the 
corresponding fitting curves for the different slab thick- 
nesses, our results are plotted in Fig. ||. The best fits 
were obtained with a bulk elastic constant Kajk^T w 66 
and an extrapolation length A/a « 2.3. The theoretical 
fitting curves agree well with the simulation results, for 
small values of k z , as one would expect for a theory valid 
for long wavelength fluctuations. At higher k z , the struc- 
ture (emphasized in the inset of Fig. |4| by a multiplying 
factor k^) is not perfectly reproduced, but the agreement 
is satisfactory. This is not surprising, since we expect 
the elastic theory to become less accurate at higher k z . 
Finally, we note that the extrapolation distance relative 
to the simulation wall position is L w + A « 2.89a, which 
compares moderately well with the value A « 2.33a ob- 
tained in the previous study of this system |28| . It should 
be noted that the director configuration of that work does 
not allow one to determine, separately, L w and A, so the 




FIG. 4: Director fluctuations (arbitrary units) as function of 
wavevector (normalized by the molecular semi-axis length a) 
for different wall separations. Symbols: Monte Carlo results; 
error estimates are indicated at some representative points, 
at higher wave-vectors the errors are smaller than the plot- 
ting symbols. Solid lines: elastic theory, fitted to parameters 
discussed in the text. Inset: fluctuations multiplied by (k z a) 2 
to emphasize structure at higher wavenumbers. Successive 
curves are offset by 0.5 for clarity. 

quoted value of A really represents A + L w . 

We have to point out some possible limitations of the 
method. The first one is computational time. The 
effect of the surfaces is to dampen the amplitude of 
long-wavelength modes; these have the longest relaxation 
times, according to eqn (|f(]), and so it is essential to carry 
out very long runs to adequately sample them. We have 
paid some attention to estimating the error bars on the 
measured values of \Q az (k z )\ 2 , as indicated in Fig|J: the 
larger values at low k z follow directly from this effect. We 
shall return to examine the time-dependence of fluctua- 
tions in a future publication. The second limitation is the 
actual sensitivity of the measured averages to the vari- 
ation in the anchoring strength and cell thickness. One 
might expect that in practice it is not possible to mea- 
sure large values of the anchoring parameter £ = L z j\ 
so we need reasonably thin cells. As the cell thickness 
L z becomes large, the fluctuation spectrum becomes in- 
sensitive to L w . However, it is important that the walls 
do not become too close: the bulk region should be suffi- 
ciently large compared to the interfacial region. Only in 
this case can we assume that the scalar order parameter 
Q in the liquid crystal bulk is constant for large scale 
fluctuation modes. 

To summarize, analysis of the director fluctuations in 
nematic liquid crystal slabs allowed us to measure the 
surface anchoring strength parameter. The method has 
been tested for a system of hard ellipsoids of revolution 
of elongation e = 15 confined between hard walls with 
homeotropic anchoring. Careful analysis of fluctuations 
in slabs of different thickness has allowed us to resolve, for 
the first time, the position of the elastic boundary condi- 
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tion relative to the simulation wall, as well as measuring 
the extrapolation length. The clastic theory gives a good 
description at low wavenumbers, but is less accurate at 
higher wavenumbers. 
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Fluctuation spectrum 

To obtain the spectrum of the q z modes of the fluctu- 
ations we need to select from the eigenfunctions of the 
operator —^KV 2 those which satisfy the boundary con- 
ditions (I). Substituting eqn (|^) into the boundary con- 
ditions (4), we obtain a set of linear equations 



(6 + ix)e ix <$n (+) + (f - ix)e" ix < 5n ( - ) 
(£ - i X )<5n (+) + (£ + i x )Jn<-> 



This set of linear homogeneous equations for 5n^> has 
a nontrivial solution if its determinant equals zero. This 
condition leads to the secular equation for the q z vector 
@ and relation (0). 

In the case of strong anchoring, £ — > oo, the q z spec- 
trum is equidistant: q z L = nn, where n is a positive 



integer. For finite anchoring coefficient £, we have a shift 
in this spectrum. The magnitude of the shift depends on 
the anchoring parameter £. Indeed, for sufficiently strong 
anchoring parameters, £ >> 1, asymptotically: 



q z L 



2nn 



1,2,. 



n/£ « 1, 



which is equivalent to replacing L by (L + 2A), A being 
the extrapolation length. 
For weak anchoring, £ << 1: 



q z L 

q z L 



2£ 



n 



1 ' 2 n = . 



1,2,. 



The spectrum of the q x ,q y wavevectors depends on the 
system geometry. Again, if we have periodic boundary 
conditions in x and y directions, the q x and q y have a 
discrete spectrum, on a fine grid q a L a = 2Tm a , (n a — 
0,1,2,...,), otherwise qj_ = {q x , qy) are unrestricted. 

It is also easy to show that the eigenfunctions which 
correspond to different eigenvalues q z and q' z are orthog- 
onal. Indeed, using eqn.(^) we can rewrite 

$(qx,?,) = <5n(+)(q ± , 9z )e i ^^+W-)(q ± , (?2 )e- i ^^ 
[X C0S (Qzz) + £ sm(q z z)} <W +) (qj_, q z 



iX + £ 



It is easy to check using the secular equation (||) that 
functions 4>(q z ) = x cos {lzz) + £ sin(q z z) are orthogonal, 
i.e. 



5 )c/>(q' z )dz = (2£ + X 2 



+x 2 +e)s, 



9*>9j 



(16) 



where Sq z ,q> is the Kronecker delta. Therefore, the eigen- 
functions $(qi, q z ) are orthogonal and can be normalized 
using eqn (fl6|). 
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